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ABSTRACT 

Stars with a different vertical motion relative to the galactic disk have a different 
average acceleration. According to Modified Newtonian Dynamics (MOND) theories 
they should therefore have a different average orbital velocity while revolving around 
the Milky Way. We show that this property can be used to constrain MOND theories 
by studying stars in the local neighborhood. With the hipparcos dataset we can 
only place marginal constraints. However, the forthcoming GAIA catalogue with its 
significantly fainter cutoff should allow placing a stringent constraint. 
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1 INTRODUCTION 

The general inconsistency between observed luminous mat¬ 
ter and the larger gravitationally inferred mass has led to 
the postulated existence of dark matter. An alternative ex¬ 
planation was proposed by Milgrom (1983), by modifying 
the Newtonian dynamics (MOND). 

According to the standard non-relativistic formulation 
of MOND (Milgrom 1983), a (Newtonian) force F acting on 
a body will cause it to accelerate according to the modi¬ 
fied second law, F = m/i (|a|/ao) a. Here p(x) is the MON- 
Dian interpolating function, which behaves as p(a; 1) « 1 
and p.{x 1) « X, while ag is a global parameter with di¬ 
mensions of acceleration. In order to explain the observed 
flat rotation curves of galaxies, ag must be of the order 
of magnitude of the galactic centripetal acceleration, i.e. 
ag ~ 10~® cm s“^. Two common choices for the interpo¬ 
lating function are the “simple” function pi{x) = x/{l +x), 
and the “standard” function I-l{x) = x j\/l + x^ . 

This simplest description (which does not obey 
the expected conservation laws) was reformulated by 
Bekenstein & Milgrom (1984) using a modified Poisson 
equation, with a non-relativistic MONDian gravitational po¬ 
tential. The resulting non-linear Poisson equation is 

V-[M(|V<?!.|/ao)V<(.] =47rGp. (1) 

An alternative non-relativistic description is that of 
Quasilinear-MOND (QuMOND, Milgrom 2010), in which 
two potentials are defined. One is an auxiliary potential 
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which solves the linear Poisson equation, and a second from 
which the force is derived. Respectively, they are 

A(()iv = drrGp, A0 = V • [p (|V0iv|/fflo) Vijijv] . (2) 

Here p(j/) is the inverted interpolating function, which is 
related to p(a;) through p(t/) = 1/pl[x{jj)], where y = xfj,{x). 
In the Newtonian regime ir{y 1) ~ 1, whereas in the deep 
MOND regime we have v{y <C 1) « For example, the 

inversion of the standard interpolating function is v{y) = 

(l/2 + {y-^ + l/Ay^^y\ 

Within this description, the MONDian gravitational 
field g = —V<?!) can be expressed in terms of the Newtonian 
field f = —V(f>N as 

g = i'(//ao)f+ cr. (3) 

This equation is derived immediately from eq. 2, where 
(7 is a curl field (V • ct = 0) that is generally non-zero and 
is defined so that V x g = 0, as required for a conservative 
(irrotational) field^. It is common to use eq. 3 while neglect¬ 
ing the curl field because this leads to a simple algebraic 
relation between the MONDian and Newtonian fields, and 
it does not involve solving a non-linear differential equation. 
Indeed, Brada & Milgrom (1995) have shown that such a 
curl field is usually small in comparison to the Newtonian 
field. In our present work we will be interested mostly in 
the vertical derivatives of the gravitational field, and it is 


^ This can be more intuitively understood by considering tr as a 
“magnetic” field generated by an effective current density jefj = 
(c/47r)(p'/ao)f X V|f|. 


© 0000 RAS 



2 B. Margalit & N. J. Shaviv 


not straightforward a priori that these derivatives are neg¬ 
ligible. For this reason, we use the full form of eq. 3 in our 
analysis. 

The non-relativistic description of MOND was also 
generalized (Bekenstein 2004) and written as a Lorentz- 
covariant theory of gravity, called the tensor-vector-scalar 
(TeVeS) theory. We will not discuss relativistic aspects, as 
they are irrelevant for the present work. 

In addition to the study of MOND in the context of 
galactic rotation curves, which were the trigger for the devel¬ 
opment of MOND, MOND theories were studied and tested 
in various additional settings. These include globular clus¬ 
ters (Scarpa et al. 2011; Ibata et al. 2011), dwarf spheroidal 
galaxies (Angus 2008), and galactic clusters, including 
both virialized and interacting (Pointecouteau & Silk 2005; 
Angus et al. 2007). Our analysis will focus on testing MOND 
using the Milky-Way rotation curve in conjunction with the 
local stellar dynamics. 

When treating stellar dynamics within the galaxy, it is 
important to remember that stars have a random compo¬ 
nent in addition to their galactic orbital motion, because of 
which they tend to oscillate radially and vertically about 
an orbiting guiding centre. Within the context of MOND, 
this additional motion offers the opportunity to locally sam¬ 
ple stars experiencing a different average acceleration. The 
typical orbital acceleration at the solar galactocentric ra¬ 
dius is of order ao. The vertical acceleration is of order 
~ 3 X 10“®(t/60 Myr)“^(2niax/100 pc) cm s“^ 
where ^max, oj and r are the amplitude, frequency, and pe¬ 
riod of the vertical oscillations. Thus, both vertical and or¬ 
bital accelerations are comparable. 

Within the simplest description of MOND (or its 
QuMOND formulation) the implications of this additional 
acceleration component are straightforward, as we will show 
in §2. Stars having a larger vertical amplitude should witness 
under MOND a smaller gravitational force, and therefore 
have a smaller average acceleration to remain in approxi¬ 
mate circular orbits. This would then generate a systematic 
drift, linking the vertical amplitude with the extent of this 
local velocity drift. We use this fact to develop a test of 
MONDian dynamics in §4. In §4.1 we discuss the hippar- 
COS astrometric data set that we use for our analysis, and in 
§4.2 search for this effect in the presently available data, but 
do not find a vertically dependent drift. This is then used to 
place actual constraints on the MOND parameter space in 
§4.3. We then discuss the limitations of this method in §5. 


2 AN AVERAGE MONDIAN “DRAG” FROM 
THE VERTIGAL ACCELERATION OF 
STARS 

The essential point of the modified dynamics that we use in 
our analysis is directly apparent in the non-potential formal¬ 
ism of MOND. Only the absolute magnitude of the acceler¬ 
ation enters as an argument of the interpolating function. 
This feature effectively generates a coupling between the 
different axes of motion of a particle, and can in principle 
be used to distinguish between Modified and Newtonian dy¬ 
namics. We begin by showing that this behaviour is a feature 
of the potential formalisms as well. Throughout the rest of 


this paper we adopt a Galactocentric cylindrical coordinate 
system (r, 9, z). 

In addition to their circular motion around the galac¬ 
tic centre, stars exhibit vertical and radial motions about 
the galactic plane and their guiding radius respectively. For 
kinematically cool disk stars, this motion is well described by 
the epicyclic approximation (see Binney & Tremaine 2008, 
chapter 3.2), in which it is assumed that the vertical and 
radial oscillations are harmonic. 

The MONDian drag effect can be qualitatively under¬ 
stood by comparing the radial (centripetal) acceleration of 
two hypothetical stars revolving in circular orbits about the 
galactic centre. The first resides in the galactic plane (2 = 0), 
while the second exhibits a vertical motion about the plane. 
Whereas the first star will experience a constant acceleration 
Or throughout its motion, the second star will experience 
a larger average net acceleration |a| = \/ap + (aI) > Ur- 
Following the MONDian argument, at smaller accelerations 
the “effective” force that a Newtonian dynamicist would in¬ 
fer, F/fi{\a\/ ao) , is increased. Therefore, the first star will 
exhibit a larger centripetal acceleration than its vertically 
oscillating counterpart. This azimuthal drift of vertically os¬ 
cillating stars relative to the non-oscillating baseline will in¬ 
crease with larger vertical accelerations a^, and hence with 
greater vertical amplitudes 2max. 

We can find an analytic estimate for the effect 
AV(zmax) using perturbation theory. Assuming a circular 
orbit for the guiding centre, and expanding the radial dis¬ 
tance and circular velocity as: R = Ro + AR, V = Vo + AV 
for a population with the same angular momentum, the two 
deviations are to first order related through 

AR = -^Ro. (4) 


We can now expand the radial acceleration in terms 
of the velocity drift AV 


(Vo + AV)" _ _Vl _ 

Ro -t- AR Ro Ro 


( 5 ) 


Next, we can also expand the gravitational field g about 
the unperturbed guiding centre at r = Ro and 2 = 0, where 
the derivates are evaluated, giving 

gr {R, z) ~ gr|o + {drgr)o \ {z^) , (6) 

where ( 2 ^) is a time average evaluated over a vertical orbital 
cycle. 

By identifying that gr|g = — Vq^ / Ro and using eq. 4, we 
find once we compare it to the MONDian acceleration that 

+ {drgr)oRo^ ^ i (d!gr), (z^) - (7) 


This equation provides the velocity drift AV of a star in 
terms of various derivatives of the MONDian gravitational 
field and the ( 2 ^) average of the star. 

We continue by deriving explicit expressions for the 
MOND field terms, which can then be plugged into eq. 7. 
Working with the QuMOND formalism, these terms can be 
found by taking derivatives of the MOND field as defined 
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by eq. 3. We make use of the following useful relations: 

dif = J ifrdifr + fzdifz) , 

/z|o ~ ~ = 0, 

/lo = -Mo- (8) 

Taking the curl of eq. 3 and using the fact that both g 
and f are irrotational fields, we find that 

drCTz — dz(Tr = — ^ [fC dz fr + fr fzdz fz 

0.0 f '■ 

-frfzdrfr - fz^Or/z] , (9) 

where C is the derivative of the interpolating function. 

Taking the vertical derivative of this equation at the 
unperturbed guiding centre, we find 

dzOr\f^ = + (drfr) {dz fz) j |/r|]^ 

+ dzdr(Jz\f^. (10) 

Using the poloidal quadrupole morphology of the curl 
field (e.g., as visualized in fig. 2 of Brada & Milgrom 1995), 
we note that dzdrOz\Q = dr{dzOz)\Q > 0. Since this term is 
positive, we can place a lower bound on dlgr\Q by discarding 
the last term. We will show below that such a bound will 
prove useful, and that in any case the leading contribution 
to the vertical derivative of the gravitational held arises from 
the algebraic part of eq. 3. This algebraic derivative is 

dlWU/oo)fr\o = vdlfr + x (11) 

Oo 

[dlfr-{dzfzf /\fr\]^. 

By combining eqs. 10 and 11, we hnd a lower bound on 
the vertical derivative of the MONDian gravitational held: 

{dlgr), > {dlfr) - ^ [{dzfzY + {drfr) {Ozfz)] . ( 12 ) 

Here v and v' are evaluated at the guiding centre. 

Plugging eq. 12 into eq. 7, we hnd an explicit solution 
to the velocity drift AV (which is negative) as a function 
of the Newtonian gravitational held f, and the interpolating 
function v. It is 

|AU| ^ v (d'^fr) - v'/oQ [{dzfz)'^ + {drfr) jdzfz)] {z^) 

Vo ^ / Ro — Rodrgr 2 

(13) 

The second and third terms in the numerator express the 
purely MONDian ehect due to the additional average ver¬ 
tical acceleration (as discussed in the beginning of this sec¬ 
tion) . On the other hand, the hrst term in the numerator and 
the second in the denominator, both of which increase the 
velocity drift, are due to the form of the galactic potential 
and they exist also in the Newtonian limit. 

Within the Newtonian terms, d^fr represents the falloff 
of the radial held perpendicular to the galactic plane. It 
physically expresses the fact that a star will spend more time 
outside the galactic plane due to its vertical oscillations, and 
thus in an average weaker radial gravitational held. This in 
turn causes a drift in the rotational velocity. 

Similarly, the term drgr describes the falloh of the held 
with increasing galactic radius. Notice that gr in this term is 


the total acceleration (or MONDian) held, and not the New¬ 
tonian gravitational held. This is useful because the term is 
kinematically well determined. Any velocity drift must be 
accompanied by an outward shift in radius (as given by eq. 
4), where the galactic potential is weaker, thus contributing 
further to the velocity drift. 

It is important to note that both these ehects are due 
merely to the galactic potential model, and exist irrespective 
of whether the motion in this potential is assumed to follow 
Modihed or Newtonian dynamics. In particular, these effects 
should contribute to give a Newtonian velocity drift, which 
is given mathematically by taking the Newtonian limit of eq. 
13, i.e., 1 / = 1,1/' — 0. As we now show, for most physically 
reasonable parameters, the Newtonian drift is signihcantly 
smaller than the expected MONDian effect. 

Eq. 13 can in principle be used to asses the velocity drift 
for any galactic mass model by simply evaluating the grav¬ 
itational field derivatives in this equation for the respective 
model. While this may prove useful for certain studies, in 
our present work we wish to estimate the velocity drift as 
broadly as possible instead of focusing on specific galactic 
models. To this extent we will try to further simplify eq. 13, 
relying on observational constraints and a few relatively well 
established assumptions regarding the local galactic held, 
predominantly the vertical harmonic and radial epicyclic ap¬ 
proximation. 

The leading term contributing to the velocity drift ex¬ 
presses the affect of an additional vertical acceleration com¬ 
ponent. Under the harmonic approximation, this term can 
be written as 

dzfz ~ -0? jv, (14) 

where cu is the vertical oscillation frequency. 

The term in the denominator, which expresses the radial 
falloff of the the galactic held, is fairly well constrained by 
observations and can be written using the Oort constants A 
and B (Feast & Whitelock 1997) as 

drgr = ^ + 2{A + B)^ « 9 X IQ-^yr-^ (15) 

This result coincides relatively well with the naive order of 
magnitude estimate drgr ~ —gr/Rd ~ where Rd is the 
disk scale length (typically in the range 2 — 4 kpc) and is 
the galactic rotation frequency. 

The vertical falloh of the galactic held (the hrst term in 
the numerator of eq. 13) is less stringently constrained by 
observations, but using the fact that the gravitational held 
is irrotational, we can roughly asses its magnitude to be 

dlfr = dr {dzfz) ~ ~ 10"® Myr“^ pc"^ (16) 

rid 

Using this order of magnitude estimate we can compare 
the contributions of the terms in the numerator of eq. 13. 
We hnd that the second, purely MONDian, term is a factor 
of about ~ up' Rd/ao > 5 larger than the hrst. This leads to a 
signihcantly more prominent ehect under Modihed dynamics 
than with Newtonian dynamics. On the other hand, the last 
term (which too is MONDian) is smaller than the leading 
(MONDian) term by roughly a factor of ~ uP jVp « 12. We 
will nevertheless retain this term in our analysis as it further 
tightens the constraints we will hnd. 

We conclude this section by giving a complete expres¬ 
sion for the velocity drift, as a function of Zmax, by averaging 
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over the vertical oscillations (which are much faster than the 
galactic rotation) and taking = z^ax/2. We also take 
into account deviations from the harmonic vertical motion 
as these will dampen the velocity drift at high «max- This 
is achieved by interpolating between harmonic acceleration 
close to the galactic plane, and constant acceleration far 
from the plane in the following manner 


=UJo /\/l + (2n 


^/Cf 


(17) 


We have here introduced a parameter with dimensions of 
length, which can be interpreted as the disk scale height. The 
Newtonian vertical field derivative can then be expressed as 
(9z/4 « - I V. 

Finally, we rewrite v and v'' in terms of jj, and ji. This 
is convenient since the argument of the latter interpolating 
functions is the well constrained galactic acceleration = 
Vq IRo, whereas the argument of the former is the unknown 
Newtonian field. Substituting these results into eq. 13, we 
arrive at a final expression for the velocity drift as a function 
of JZmax. It is 


|AF| > 
Vo 


(dlfr) + 


Wq/aq 


p + {V^/Roao)fi' 1 + (2n,ax/C)^ 


+ 


2 

UJo 






■ 8V^/Ro - SVh (A + B)’ 


(18) 


In the following section we compare this analytic for¬ 
mula (or rather a variant of it excluding the curl field) 
against the results of numeric simulations and show that the 
results are in good agreement. We later use this expression 
in §4.3 to constrain the parameters so that AV is consistent 
with observations. 

In developing our results we have assumed that the stel¬ 
lar population is kinematically cool such that it is well de¬ 
scribed by the radial epicyclic approximation. This is an im¬ 
portant point, as our theoretical prediction pertains to these 
young disk stars alone. Halo stars are not described by our 
model, but should “seem” to drift relative to the Local Stan¬ 
dard of Rest (LSR) if included in a sample of analyzed stars. 
This extraneous effect would arise since halo stars have a 
much smaller (z-component) angular momentum than disk 
stars, and would appear to have a V component drift ve¬ 
locity relative to the disk stars. It is therefore essential to 
remove any halo stars from the data set before searching for 
any velocity drifts. This is discussed further in §4.1. 


3 THE DRAG IN A REALISTIC POTENTIAL 

Eq. 18 was obtained under a simplified potential. In this sec¬ 
tion we numerically estimate the expected drag for a realistic 
Galactic potential, as well as verify our analytic results. 

We simulate test particle motion in a static Galactic 
potential under both Newtonian and MONDian dynamics. 
We first simulate the motion under a specific, yet complete. 
Galactic mass model, so as to asses the general expected 
magnitude of the drift effect. We model the Galactic poten¬ 
tial in this case by integrating the density distribution given 
by ‘Modell’ in Dehnen & Binney (1998). For the MOND 


simulation we take only the baryonic contribution to the 
density distribution, while for the Newtonian simulation we 
include also the dark matter halo as described by the full 
model of Dehnen & Binney (1998). To obtain the MOND 
potential we use the purely algebraic relation of eq. 3 (while 
neglecting the a curl field for computational simplicity). 

We use the standard interpolating function, and obtain 
a self consistent value for aq by “fitting” the MONDian rota¬ 
tion curve to its Newtonian counterpart. This method yields 
a value of aq = 2.06 x 10~® cm s“^ for the MOND acceler¬ 
ation parameter. Although this value is slightly higher than 
the 1.2 X 10“® cm s“^ found from surveys of external galax¬ 
ies (Begeman et al. 1991; Gentile et al. 2011), it has already 
been pointed out that higher values of ao are required to fit 
the Milky Way rotation curve using the standard interpo¬ 
lating function (e.g., Famaey & Binney 2005). 

We then integrate the equations of motion for a series 
of test particles having the same initial Galactic radius Ro 
and velocity Vo, but with different values of Zmax. This is 
achieved by implementing a 2D leapfrog integration scheme 
in the R-z plane. The radial field component is modified to 
account for rotation by adding the effective term L^/R^, 
where Lz is the test particle’s specific angular momentum. 
For each run of the simulation, we extract the guiding ra¬ 
dius velocity (V) about which the particle oscillates verti¬ 
cally and radially. The velocity drift AH is then found by 
subtracting Vo from each of these (V) . The results compar¬ 
ing the Newtonian and MONDian velocity drifts obtained 
by these simulations are shown in figure 1 (b). It is clear 
from the figure that as expected in both cases the lag in¬ 
creases with larger 2max, and that the effect is substantially 
more dramatic in the MOND case. Additionally, while the 
Newtonian lag is very nearly quadratic in Zmax, the MON¬ 
Dian effect breaks from this trend, due to non-Gaussianity 
in the vertical potential. 

We also consider the simplified model discussed in ob¬ 
taining Eq. 18. This allows us to verify the validity of our an¬ 
alytic results, as well as present the dependence on the var¬ 
ious parameters. For the vertical field component we take a 
linear model at small ^max, with the slope Wq, which evolves 
into a constant field at larger Zmax, exactly as presented in 
the previous section. For the radial component we take a 
linear dependence on radius and quadratic dependence on 
vertical displacement, giving us overall: 


1 + (^max/ C) 

fr = fro + {drfr)o {R - Ro) + ^ 

We once again simulate MONDian and Newtonian test 
particle motion in this field and extract AH(2inax)- We plot 
the numerical results as well as the analytic expression for 
the velocity drift (a revision of eq. 18 excluding the a field 
contributions). The simulation results agree with the theo¬ 
retical expression to within a few percent, with larger devi¬ 
ations from the analytic approximation at larger ^max- The 
main source of deviation from eq. 18 is terms of higher or¬ 
der in (a^)/a^ which were neglected in deriving the result 
for sake of simplicity and coherence. Figure 1 (a) shows ex¬ 
amples of the simulation results for a few different param¬ 
eter values, as well as the respective analytic curves. One 
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Figure 1. The V velocity drift as a function of the vertical oscillation amplitude ^max- (a) Simulation results (different symbols) 
and analytic approximation (solid lines) for the “simplistic” galactic model described in the text. The analytic approximation is in 
good agreement with the simulation results, and begins to diverge slightly only at larger ^max as expected. The MONDian effect is 
characteristically several times larger than the expected Newtonian drift, but can be quenched at higher 2:max by taking a small non¬ 
linearity scale g (b) Simulation results for the galactic potential obtained by integrating mass model ‘Modell’ (from Dehnen & Binney 
1998). The results are similar in nature to those of panel (a). The velocity drift breaks from its quadratic dependence on 2max at around 
~ 100 pc, signifying a deviation from a purely harmonic vertical potential. 


can notice from this figure the affect of choosing smaller ( 
values—effectively quenching the additional vertical accel¬ 
eration component, and hence the velocity drift, at higher 
^max- An additional way to decrease the velocity drift would 
be choosing very small values for the MOND acceleration 
parameter ao, as in this limit we revert back to classical 
Newtonian dynamics. 


4 A TEST FOR MONDIAN DYNAMICS 

In previous sections we have presented the physical argu¬ 
ments and theoretical predictions of a MONDian velocity 
drift. In the following we propose using available astrometric 
data to plot the rotation velocity Y as a function of vertical 
oscillation amplitude 2max for a selection of stars in the solar 
neighborhood. By checking for a systematic lag AY(2inax) 
we attempt to constrain parameter values under MONDian 
dynamics. 

We assume a harmonic vertical potential and assign a 
vertical amplitude to each star in our sample using the star’s 
observables W (vertical velocity) and z (vertical displace¬ 
ment), that is 

2max = (20) 

The harmonic assumption is valid as long as we con¬ 
strain ourselves to small vertical amplitudes, we therefore 
include anharmonic effects in our model as discussed in 
the previous section. Bahcall & Bahcall (1985) show that 
the expected non-Gaussianity of the potential is ~ 15% at 

Z^max = 100 pc. 


4.1 The astrometric data 

We use the Extended Hipparcos Compilation (XHIP) cata¬ 
log (Anderson & Francis 2012) in our analysis. This recent 
catalog is based on the New Reduction of Hipparcos data 
(van Leeuwen 2008) which includes parallax and proper mo¬ 
tion entries for 117,955 stars in the solar neighborhood. The 
catalog also uses proper motions from the Tycho-2 Cata¬ 
log (H0g et al. 2000) in cases where the Hipparcos proper 
motions exceed Tycho error bounds by some extent. Addi¬ 
tional, the catalog assigns radial velocities to 46,392 stars in 
the sample by cross referencing 47 sources of radial velocity 
data (e.g. Gontcharov 2006). We restrict our sample to kine¬ 
matically complete well localized stars by taking only stars 
for which both proper motions and radial velocity exist, and 
by imposing a cutoff for stars with a quoted parallax error 
greater than 20%. Thus we are left with a total of 32,418 
stars in our sample. 

In order to avoid possible contamination by halo stars 
we exclude relatively old stars for which B — V < 0.68mag. 
We additionally exclude high velocity outliers by restricting 
our sample to stars within the velocity ellipsoid (17/200)^ -I- 
(Y/80)^ -I- (lT/120)^ < 1. The constraints on U and W 
are quite loose for this cutoff (au ~ 25 km s“^ and 
aw ~ 11 km s“^), yet the limit on W is essentially im¬ 
material as our method explicitly uses only specific ranges 
of W (in the process of calculating Zmax) which are well 
below this limit, and the cutoff on U was found to have lit¬ 
tle significance on the final results. The V cutoff on the 
other hand is more stringent but is still well justifiable, 
since ay ~ 13—17 km s“^. Throughout the analysis we 
take the LSR as (H, Y, 1T)lsr = (7.5,13.5,6.8) km s”*^ 
(Francis & Anderson 2009). 
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4.2 Application of the Method to the Extended 
Hipparcos dataset 

We focus now on applying our method to the sample dis¬ 
cussed above. A non-trivial task in performing this anal¬ 
ysis is applying the statistics to an unbiased population. 
Since the number density of stars increases with decreasing 
galactic radii, we generally expect to find a non-Gaussian 
distribution of V exhibiting an asymmetric drift. This arti¬ 
fact can be remedied by using only stars with a particular 
galactic angular momentum L^. Since a given value of Lz 
corresponds to a particular guiding radius (for a population 
with some 2max), stars with the same angular momentum 
will be observed in different phases of epicyclic oscillation 
about a nearly fixed guiding radius, and should be sampled 
from the same number density at this galactic radius. 

For this purpose we restrict our sample to stars with 
galactic angular momentum in the range of ALz = 1.5 x 
10"* pc km s“^ around Lz ~ 1.85 x 10® pc km s“^. This value 
was chosen to avoid known streams, and samples 1109 stars 
in the solar neighborhood with (V) ~ Vq. The remaining 
kinematic distributions are therefore completely Gaussian 
and unaffected by any streams or an asymmetric drift. The 
range of ALz was chosen such as to maximize the number of 
stars in the sample without affecting its inherent standard 
deviation (which is on the order of av ~ 2.5 km s“^), in 
other words, ALz < Rav- We note that even without this 
extra precaution we obtain virtually the same constraints 
on MOND. This implies that the possible systematic errors 
from a biased population are not important in the present 
dataset, and that the reduced noise in this subset roughly 
cancels the smaller sample size. Though, it should be men¬ 
tioned that the subset appears to have more Gaussian dis¬ 
tributions. 

We first calculate jZmax for each star in our sample 
using equation (20) under the mentioned assumption of 
To — 64 Myr. We divide the sample obeying ^max < 256 pc 
into 5 consecutive Zmax bins containing an equal number 
of stars (209 each). For each bin we then calculated the 
mean rotational velocity V by fitting the Gaussian distri¬ 
bution, and plotted the results as a function of Zmax (taken 
as the mean 2max for the population of each bin). Figure 2 
shows the resulting plot of Al/(2max) (normalized such that 
AV = 0 for the first bin), as well as the underlying veloc¬ 
ity histograms for each bin. Since the distributions are very 
nearly Gaussian, using the arithmetic mean instead of the 
guassian fit results in nearly identical results. 

The figure also depicts the drift velocity for the sub¬ 
population of stars moving away and towards the galactic 
plane. A significant discrepancy between the two could arise 
if the sample of stars have a non-negligible contamination 
of stars which haven’t been equilibrated, or if the sample is 
biased by the inclusion of a prominent stellar stream / mov¬ 
ing group. The fact that no significant difference is evident 
between the two subsamples indicates that our population 
is not affected by such contamination. 

When older stellar populations were included in the 
analysis by removing the color cuts and when the restriction 
on Lz was lifted, the V distributions were far from Gaussian 
and showed a prominent asymmetric drift. In this case, the 
AV (zmax) (calculated this time as the arithmetic mean of 
each bin) exhibited a negative trend of order ~ 2 km s“^ 


at 2inax ~ 200 pc. This behavior is as expected from a con¬ 
tamination of the sample by an asymmetric drift, which is 
accentuated once older stars are introduced into the sam¬ 
ple (older stellar populations exhibit larger epicyclic motions 
and can originate from smaller galactic guiding radii). Older 
halo and thick disk stars, which may or may not be bound 
to the galactic plane and are therefore not well described 
by our model can also contaminate the sample in this case. 
Although halo stars may only amount to ~ 0.15% of the 
thin disk population, stars associated with the thick disk 
can contaminate our sample more significantly. As shown 
by Girard et al. (2006); Moni Bidin et al. (2012), these stars 
exhibit a substantial shear of 20 — 30 km s“^ at the galac¬ 
tic plane, and their number density is a non-negligible 10% 
of thin disk stars. Indeed, we find a clear indication of this 
thick disk population at V ~ —25 km s“^ when the Lz cut 
is removed. 

4.3 Constraining MOND parameters with the 
observational results 

Using the observational results obtained in the previous 
section along with the theoretical expectation for a MON- 
Dian drift, we can constrain the MOND parameter space. 
The data analysis as summarized in figure 2 is consis¬ 
tent with a flat AV curve. As such, and because the 
velocity drift is monotonic in 2max, we can restrict any 
velocity drift at 2max = 190.7 pc to be no more than 
0.375 km s“^ (0.823 km s“^) at Icr {2a) confidence levels re¬ 
spectively (which is the value of the furthest point in figure 
2 including uncertainties for the initial, “zero” 2max point). 

In the following analysis we make use of eq. 18, tak¬ 
ing the values of the Oort constants as A = 14.82 ± 
0.84 km s“^ kpc“^,i3 = —12.37 ± 0.64 km s“^ kpc“^ 
(Feast & Whitelock 1997). We also considered that the ob¬ 
served vertical oscillation frequency assuming Newtonian 
dynamics is uo = 27r/64 Myr“^ (Shaviv et al. 2014). Note 
that a longer period would scale down the vertical accelera¬ 
tions and with it the constraints placed on MOND, inversely 
proportional to the fourth power of the period. 

Additionally, for the solar radius we take Ro — 8 kpc 
(Reid 1993). The galactic rotation velocity is then given as 
Vo = Ro{A — B) = 217.52 km s“^. For the vertical falloff 
of the radial field a typical order of magnitude value may 
be used, as given by eq. 16. This value is consistent with 
the value obtained by analysis of the galactic mass model of 
Dehnen & Binney (1998) ‘Modell’ (this model is discussed 
more in §3). We are are left with two parameters we wish to 
constrain, oq and ({, as well as the choice of the interpolating 
function when evaluating eq. 18 at a given 2max. The veloc¬ 
ity drift can always be decreased by taking smaller values 
for the MOND acceleration ao, because in the limit where 
ao —S' 0 we revert back to a Newtonian theory. Similarly the 
velocity drift can be decreased by taking a smaller nonlin¬ 
earity scale height C,, effectively reducing the average vertical 
acceleration. 

We thus leave ao and C, as free parameters, and for a 
given interpolating function ^{x) relate the two by solving 
numerically the constraint: 

AU [ao, C] (190.7pc) < 0.073 ± 0.448 km s“^ (21) 

We solve eq. 21 for the two commonly used interpolat- 
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0 < Zjjiax ^ 50 50 < Zjiiax 77.7 77.7 < Z^ax < 110 110 < Zjnax 150 150 < Ziiiax 256 




Figure 2. The average V-component deviation of the selected stars as a function of the amplitude of their vertical motion (with respect 
to the galactic plane). Panel (a) shows the V distributions for the 5 Zmax bins. Panel (b) shows the resulting deviation of each bin’s 
mean V from the first bin. The resulting curve is centered around zero, and shows no appreciable trend. Specifically, no lag is evident. 
The black errorbars are the statistical standard error of the mean calculated for each bin. The upwards (and downwards) facing triangles 
denote the drift velocity of the subpopulation of stars moving away (towards) the plane. 



Figure 3. Constraints on the (ao,C) parameter space for any MONDian theory. The constraint is found by evaluating eq. 18 and 
requiring that AV^ [ao? C] (190.7 pc) < 0.073di0.448 km s~^ to be consistent with observations (as presented in §4.2). Figure 3(a) displays 
the evaluated parameter space for the ‘simple’ interpolating function, whereas figure 3(b) displays results for the ‘standard’ interpolating 
function. The regions to the right of the black solid/dashed curve in both figures illustrate the excluded regions to within one/two standard 
deviations (respectively), and taking a conservative lower bound for the vertical falloff of the radial gravitational field, d^fr = 0. Similarly, 
the regions to the right (and above) the red solid/dashed curves show the 1/2(t parameter space additionally excluded by taking a more 
realistic constraint of d^fr = 10“® Myr“^ pc“^. 


© 0000 RAS, MNRAS 000, 000-000 




















8 B. Margalit & N. J. Shaviv 


ing functions, in each case giving both a loose constraint by 
taking d^fr = 0, and a more stringent constraint by choos¬ 
ing a realistic value of d^fr ~ 10~® Myr~^ pc“^. The results 
are plotted in Fig. 3 for the standard and simple interpolat¬ 
ing functions. The regions to the right of the black solid 
(dashed) curves in both figures illustrate the excluded re¬ 
gions to within one (two) standard deviations (respectively), 
for the loose constraint (i.e. taking d^fr = 0). Similarly, the 
regions to the upper-right of the red solid (dashed) curves 
show the la {2a) parameter space excluded by taking the 
more realistic estimate for fr¬ 
it is worth noting at this point that the parameter ({ 
is essentially the vertical scale height for the galactic disk. 
Traditional mass models (for e.g. Dehnen & Binney 1998; 
Bissantz & Gerhard 2002) take this parameter as about 
185 pc for the thin disk (with which we are dealing). Ad¬ 
ditionally, we can use the observation of Bahcall & Bahcall 
(1985) that the harmonic period of vertical oscillations dif¬ 
fers by about 15% (possibly less) from the actual period at 
2max = 100 pc. lu our case the actual oscillation period is 
dependent on the scale height ^ through eq. 17. By plug¬ 
ging in the constraint given by Bahcall & Bahcall (1985) we 
conclude once again that a consistent value for ^ should be 
about 176 pc (or more). 

From figures 3(a),3(b), we notice that for the accepted 
value of the MONDian acceleration ao = 1.2 x 10~® cm s“^ 
(Begeman et al. 1991), the conservative lower bound re¬ 
stricts g ^ 164 pc at Ict for the ‘simple’ function, a value 
which is marginally inconsistent. The more stringent con¬ 
straint on the other hand yields ^ 57 pc within one stan¬ 
dard deviation. The conservative constraint on the ‘stan¬ 
dard’ function is similar yet less limiting, because the asymp¬ 
totic y-axis value of the solid curve lies higher, and gives 
(( ^ 213 pc for ao — 1.2 x 10“® cm s“^, which is marginally 
consistent. Still, if we take into account that larger values 
of ao (as much as 2 — 3 x 10“® cm s“^, Famaey & Binney 
(2005)) are needed to fit the Milky Way rotation curve with 
this interpolating function, the conservative constraint will 
permit values of ^ < 130pc, which are below the expected 
value of the disk scale-height. The tighter constraint once 
again restricts ( to values that are inconsistent with the la 
constraint and may plausibly be ruled out. 

Also note that we have not included the affect of ( in 
the analysis of the Hipparcos data (described in § 4.2), where 
we calculated ^max under the assumption of a perfectly har¬ 
monic potential. Qualitatively, introducing the non-linearity 
in would increase all calculated 2max in fig. 2 and lead to 
tighter constraints. For our specific choice of interpolation 
between linear and constant field (eq. 17), this increase in 
2max will uot Compensate for the reduction of the average 
vertical acceleration. Therefore in the limit <(■—>■ 0 MOND 
can always be “saved” (assuming d^fr = 0). Nonetheless, 
we would expect a slight shift of all curves in fig. 3 to the 
left. In this sense we have underestimated our constraints 
on the MOND parameter space. 

Finally, we show in figure 4 the change of the 2a real¬ 
istic curves in the (ao,C) plane under variation of the as¬ 
sumed observational parameters Ro, A, B,iJo- This provides 
a measure for considering the effect of the uncertainties as¬ 
sociated with these parameters. For the Oort constants, we 
take the stated uncertainties of Feast & Whitelock (1997). 
Based on different studies, the solar galactocentric radius 


tends to vary around 8 kpc, but is almost certainly in the 
range of 7.5 — 8.5 kpc. We therefore take the variation in 
Ro as ±500 pc. Finally, we allowed the vertical oscillation 
period to vary signihcantly, by ±8 Myr. This parameter is 
probably the most important for the MONDian drag effect, 
but there is also less of a consensus regarding its value which 
is why we permit such large fluctuations. 

Each of the curves in figure 4 was obtained by chang¬ 
ing one of the parameters, while keeping the others fixed at 
their nominal values. It is clear that the Oort constants and 
galactocentric radius only modestly alter the nominal curve, 
so that the analysis is relatively robust for these parameters. 
As expected, the vertical oscillation period ro has a larger 
affect on the results. Even so, the impact of decreasing the 
period is less substantial than one would naively expect from 
plugging in larger cuq in Eq. 18. This is due to a proportional 
increase in all «max values obtained by Eq. 20, and in par¬ 
ticular that of the furthest bin in figure 2. We thus obtain 
tighter observational constraints which somewhat compen¬ 
sate for the larger uq. 

Since we “rule out” MOND at only < 2a (for the reason¬ 
able galactic potential), the analysis with the present dataset 
can only be considered as an indication, and a proof of con¬ 
cept for the method. 


5 DISCUSSION 

We have shown both theoretically and numerically that a 
rotational velocity drift should be associated with stars os¬ 
cillating perpendicular to the galactic plane, and should in¬ 
crease in magnitude with the maximal vertical extent of the 
star’s motion. This effect is inherently different under the as¬ 
sumption of MONDian versus Newtonian dynamics, and we 
have shown that the MONDian velocity drift is significantly 
larger than its Newtonian counterpart. Although our results 
were obtained under the QuMOND formulation of the mod¬ 
ified dynamics, alternative MONDian theories differ from 
this formulation (in the non-relativistic limit) only by a curl 
held such as a. Since we have found that the main contri¬ 
bution to the MONDian drift arises from the algebraic part 
of relation 3, we expect that alternative formulations will 
yield very similar results. In particular, eq. 18 excluding its 
smaller curl held contributions can be obtained exactly and 
directly from the non-linear Poisson formulation of MOND. 

Previous works by Stubbs & Garg (2005) and 
Nipoti et al. (2007) have been motivated by the same 
concept of interplay between vertical and radial dynamics 
in MOND, yet differ from our current work in several 
aspects. Specihcally, neither of these works manage to 
quantitatively constrain MOND with current observational 
data. Stubbs & Garg (2005) method is also somewhat 
limited in its applicability, as it specihcally assumes an 
exponential disk model with a non-varying scale height, 
and relies on data measured over a large range of galactic 
radii (throughout which this assumption most likely does 
not hold). Our present method on the other hand, relies 
only on local kinematic observables which can be better 
constrained. Being also well aware of the limitations of 
the harmonic approximation, we explicitly introduce and 
allow deviation from this model by proper choice of the free 
parameter which we carry throughout. 
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(b) 


Figure 4. Change of the “realistic” 2a curves in the (ao,C) parameter space under variation of the assumed physical variables within 
their plausible uncertainties. While the Oort constants A and B and the solar galactocentric radius Rg only modestly alter the nominal 
curve, the larger deviation in ro affects the result more significantly. 


By using complete kinematic data for stars in the lo¬ 
cal galactic neighborhood obtained from the XHIP cata¬ 
log, we have shown that observations are consistent with 
no apparent velocity drift up to Zmax = 190 pc and within 
our confidence bounds of ±0.448 km s“^. This fact allows 
us to constrain the (aoiC): ^ind u(^) parameter space by 
requiring that they will be consistent with this observa¬ 
tion. We use our analytic expression (lower-bound) for the 
velocity drift (eq. 18) as a function of these parameters 
and conclude that current accepted values for ag and the 
disk scale-height are marginally inconsistent with our pre¬ 
sented observational constraint. By assuming a more real¬ 
istic value for the vertical falloff of the radial gravitational 
field {d^fr = 10~®Myr~^pc“^) we find that reasonable val¬ 
ues for these two parameters can be ruled out at 68% con¬ 
fidence bounds, and conclude that present MOND theories 
may suffer from inconsistencies. Since this initial analysis 
yields only weakly significant constraints, it serves primarily 
as an indication of such possible MONDian inconsistencies, 
and should more importantly be interpreted as a proof of 
concept for our method. 

We expect that even the conservative constraint (ne¬ 
glecting d^fr) may be improved in the near future with the 
GAIA mission (Perryman et al. 2001), by minimizing the 
error bounds on AV and possibly increasing the maximal 
observed Zmax. Even a decrease of these error bounds by a 
mere factor of two will significantly impact the results as 
the theory is highly non-linear. By increasing the sample 
size and detection horizon, we may also be able to constrain 
the velocity drift at larger ^max than present, which will pos¬ 
sibly constrain the (ao,C) parameter space even further. It 
is important to note that our method was developed as a se¬ 
ries expansion about the galactic plane, and may thus break 
down at very large vertical extents. Nonetheless, even a mod¬ 
est increase in 2max to ~ 300 pc (which is likely still well 
described by our method, especially since we take into ac¬ 
count the non-linearity of the potential at larger heights via 
the parameter () would tighten the constraints on MOND 
significantly. 


An interesting by-product result of our analysis is an 
upper bound on d^fr assuming normal Newtonian dynam¬ 
ics. This bound is easily derived by taking the Newto¬ 
nian limit of eq. 21, i.e. /r = 1,/r' = 0, and solving 
for d^fr- Using this procedure we find: < 2.1 x 

10“® Myr“^ pc~^ (4.7 x 10“® Myr“^ pc~^) at Ict (2cr) con¬ 
fidence level respectively. This result may be useful in con¬ 
straining Newtonian galactic mass models, but lies outside 
the scope of our present work. 
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